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ABSTRACT 

A  method  to  simulate  digital  human  running  using  an 
optimization-based  approach  is  presented.  The  digital 
human  is  considered  as  a  mechanical  system  that 
includes  link  lengths,  mass  moments  of  inertia,  joint 
torques,  and  external  forces.  The  problem  is  formulated 
as  an  optimization  problem  to  determine  the  joint  angle 
profiles.  The  kinematics  analysis  of  the  model  is  carried 
out  using  the  Denavit-Hartenberg  method.  The  B-spline 
approximation  is  used  for  discretization  of  the  joint  angle 
profiles,  and  the  recursive  formulation  is  used  for  the 
dynamic  equilibrium  analysis.  The  equations  of  motion 
thus  obtained  are  treated  as  equality  constraints  in  the 
optimization  process.  With  this  formulation,  a  method  for 
the  integration  of  constrained  equations  of  motion  is  not 
required.  This  is  a  unique  feature  of  the  present 
formulation  and  has  advantages  for  the  numerical 
solution  process.  The  formulation  also  offers 
considerable  flexibility  for  simulating  different  running 
conditions  quite  routinely.  The  zero  moment  point  (ZMP) 
constraint  during  the  foot  support  phase  is  imposed  in 
the  optimization  problem.  The  proposed  approach  works 
quite  well,  and  several  realistic  simulations  of  human 
running  are  generated. 

INTRODUCTION 

The  3D  digital  human  running  problem  is  formulated  as 
an  optimization-based  predictive  dynamics  problem.  It  is 
noted  that  the  expression  “predictive  dynamics”  has 
been  coined  to  characterize  a  class  of  physics-based 
problems  that  are  modeled  using  differential  equations  of 
motion  and  that  would  otherwise  require  the  solution  of 
such  problems  using  time-step  integration.  However,  in 
this  case,  predictive  dynamics  is  a  broadly  applicable 
formulation  for  addressing  such  problems  using 
optimization  techniques  without  having  to  integrate  the 
equations  of  motion.  Indeed,  the  formulation  does  not 
have  a  unique  solution,  and  constraints  play  an 
important  role  in  the  solution  process. 


In  the  problem  of  running,  the  objective  is  to  predict  (or 
calculate)  joint  angles  and  torques  at  the  joints  over 
time,  also  called  joint  and  torque  profiles,  respectively. 
For  the  problem  of  running,  a  minimal  set  of  constraints 
is  imposed  in  the  formulation  of  the  problem  to  simulate 
natural  running  of  the  digital  human.  The  human  running 
problem  is  distinguished  from  the  walking  problem  in  that 
there  is  a  flight  phase  during  each  step  of  running.  In  the 
present  formulation,  running  steps  are  assumed  to  be 
periodic  and  symmetric  for  the  right  and  left  steps.  Both 
the  support  phase  and  the  flight  phase  are  modeled.  A 
companion  paper  by  Xiang,  et  al  (2007)  presents  the 
formulation  for  the  human  walking  problem. 

The  digital  human  is  modeled  as  a  mechanical  system 
that  includes  link  lengths,  mass  moments  of  inertia,  joint 
torques,  and  external  forces.  The  entire  model  has  55 
degrees  of  freedom  (DOF) — 6  for  global  translation  and 
rotation  and  49  for  the  body.  A  DOF  in  this  case 
characterizes  a  kinematics  jointed  pair  in  the  kinematics 
sense,  where  various  segments  of  the  body  are 
assumed  to  be  connected  by  revolute  joints.  The  B- 
spline  interpolation  is  used  for  time  discretization,  and 
the  Denavit-Hartenberg  (DH)  method  is  used  for 
kinematics  analysis.  The  recursive  Lagrangian 
formulation  is  used  for  the  equations  of  motion;  it  was 
chosen  because  it  is  considered  to  be  quite  efficient.  The 
equations  of  motion  are  verified  by  the  forward  dynamics 
process  using  a  commercial  general-purpose  multi-body 
dynamics  software  code. 

The  problem  is  formulated  as  a  nonlinear  optimization 
problem.  A  unique  feature  of  the  formulation  is  that  the 
equations  of  motion  are  not  integrated  explicitly;  this  has 
become  the  most  important  contribution  because  it 
provides  for  a  generalized  method  to  solve  dynamic 
indeterminate  problems  that  would  otherwise  require 
computationally  intensive  integration  methods.  They  are 
imposed  as  equality  constraints  in  the  optimization 
process.  An  algorithm  based  on  the  sequential  quadratic 
programming  approach  is  used  to  solve  the  nonlinear 
optimization  problem.  The  control  points  for  the  joint 
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angle  profiles  are  treated  as  design  variables;  when 
calculated,  they  provide  for  the  complete  motion.  For  the 
performance  measure  in  the  optimization  problem,  the 
mechanical  energy  that  is  represented  as  the  integral  of 
the  sum  of  the  squares  of  all  the  joint  torques  is 
minimized.  The  dynamic  stability  is  achieved  by 
satisfying  the  zero  moment  point  (ZMP)  constraint  in  the 
support  phase.  The  zero  yawing  moment  (ZYM) 
constraint  is  imposed  so  that  the  upper-body  motion  is 
compensated  by  the  lower-body  motion.  The  solution  is 
simulated  in  the  Santos™  human  modeling  and 
simulation  environment.  The  simulations  show  a  very 
natural  running  motion  of  the  digital  human.  The  joint 
torque  and  angle  profiles  and  ground  reaction  forces  are 
recovered  from  the  simulation.  A  couple  of  simulations  of 
running  at  different  speeds  and  carrying  loads  on  the 
back  are  generated  quite  routinely  and  efficiently. 

Virtual  human  dynamics  simulation  is  a  very  active  area 
of  research.  In  recent  years,  many  papers  have  been 
published  on  biped  digital  human  motion.  For  the  digital 
human  running  problem,  however,  there  are  only  a  few 
papers  in  the  robotics  area.  In  the  biomechanics  area, 
reports  have  been  mostly  on  experimental  research 
involving  subjects  rather  than  the  dynamics  simulation  of 
the  problem. 

Biomechanics'.  Cavanagh  and  Lafortune  (1980)  studied 
ground  reaction  forces  and  center  of  pressure  (COP) 
patterns  for  18  subjects  during  running.  Simpson  and 
Bates  (1990)  experimented  with  the  effect  of  running 
speed  on  joint  moments  in  the  support  phase.  Ounpuu 
(1994)  described  the  terminology,  kinematics,  and 
kinetics  of  human  walking  and  running  in  clinics  in  a 
sports  medicine  journal.  Novacheck  (1998)  summarized 
the  biomechanics  of  human  running.  He  discussed 
kinematics,  kinetics,  COP,  and  muscle  activities  in 
relation  to  EMG  results  and  injuries. 

Robotics :  Many  researchers  have  worked  on  the 
problem  of  walking  robots.  Since  walking  and  running 
are  considered  part  of  biped  locomotion,  there  are  many 
papers  in  the  robotics  walking  area  that  are  relevant  for 
the  running  problem.  For  example,  the  zero  moment 
point  (ZMP)  dynamics  stability  constraint  (Vukobratovic 
et  al.,  1990)  can  be  included  in  the  running  problem  as 
well  as  the  walking  problem.  However,  only  the  literature 
about  the  running  problem  will  be  discussed  in  this 
section.  Honda  has  been  developing  humanoid  robots 
since  1986,  and  their  robot  Advanced  Step  in  Innovative 
Mobility  (ASIMO)  is  the  most  advanced  running  robot  to 
date.  ASIMO  can  run  at  6  km/h.  Fujimoto  (2004)  used  an 
optimization  technique  for  creating  trajectory  for  a  biped 
running  robot.  The  idea  was  to  minimize  energy 
consumption  due  to  the  biped  robot’s  running  motion  by 
determining  the  joint  angles  and  torques.  This  work  was 
done  for  a  2D  model  that  had  only  7  DOF.  Nagasaki  et 
al.  (2003)  generated  the  running  pattern  by  using  the 
angular  momentum  and  the  control  theory;  however, 
only  a  few  seconds  of  simulation  was  performed. 
Roussel  et  al.  (1998)  published  about  the  generation  of 
an  energy-optimal  complete  gait  for  biped  robots.  This 


work  did  not  discuss  the  running  problem,  but  it  included 
an  impulse  term  in  the  cost  function,  which  will  be 
considered  in  the  running  problem  in  the  present  study. 
Park  and  Kwon  (2003)  developed  a  biped  robot’s 
running  motion  by  using  the  impedance  control.  Hybrid 
Zero  Dynamics  (HZD)  was  presented  by  Westervelt  and 
Grizzle  (2003).  The  robot  walked  with  quite  a  natural 
motion  with  that  method,  but  it  did  not  have  3D  stability. 

Computer  Graphics :  Basically,  the  methods  to  generate 
locomotion  in  robotics  and  animation  are  similar. 
However,  animators  are  more  interested  in  high-level 
behavior,  while  researchers  in  robotics  are  interested  in 
joint  torques  and  forces.  Hodgins  (1996)  simulated  3D 
digital  human  running.  The  approach  basically  used  the 
control  theory  for  the  mechanical  system,  and 
commercial  software  was  used  for  the  dynamics  solution 
of  the  mechanical  system.  However,  3D  stability  was  not 
considered  in  her  work.  Kang  et  al.  (1999)  proposed  a 
model  based  on  a  one-legged  planar  hopper  with  a  self¬ 
balancing  mechanism  for  human  running  animation. 

In  previous  work,  Kim  et  al.  (2004)  developed  an 
optimization-based  dynamic  motion  prediction  method 
for  the  digital  human  gait.  They  used  B-spline,  the  DH 
method,  and  ZMP/ZYM  stability  for  the  walking  problem. 
However,  the  joint  torques  could  not  be  obtained  with 
that  formulation  because  the  equations  of  motion  were 
not  considered.  Currently,  the  equations  of  motion  and 
their  sensitivities  are  included  so  that  the  joint  torques 
can  be  obtained. 

Considering  the  human  running  gait  cycle,  the  running 
style  depends  on  the  speed  of  running.  For  example,  at 
slower  running  speeds,  the  heel  touches  the  ground  first. 
In  fast  running  or  sprinting,  the  fore-foot  touches  the 
ground  first.  Moreover,  the  upper  body  motion  is  different 
for  slower  and  faster  running  (sprinting).  The  faster  the 
runner’s  speed,  the  more  arm  swinging  motion  is 
generated  to  minimize  energy  consumption.  Running  is 
differentiated  from  walking  not  by  the  speed  but  by  the 
existence  of  a  flight  phase.  During  a  walk,  whether  slow 
or  fast,  there  exists  a  double  support  phase  (where  both 
feet  are  on  the  ground).  The  period  from  the  initial 
contact  of  one  foot  to  the  following  contact  of  the  same 
foot  is  called  the  gait  cycle.  One  gait  cycle  of  running  is 
composed  of  two  phases:  the  support  phase  and  the 
flight  phase.  The  flight  phase  starts  with  a  toe  off  and 
ends  with  the  strike  of  the  other  foot.  The  support  phase 
starts  with  a  foot  strike  and  ends  with  same  foot’s  toe  off 
(Figure  1).  In  the  area  of  biomechanics,  the  distance 
from  one  foot’s  strike  to  the  other  foot’s  strike  is  called  a 
step.  Also,  the  distance  from  one  foot’s  strike  to  the 
same  foot’s  subsequent  strike  is  called  a  stride. 


/„,}  and  control  points  q0  qx  q2  ,  the  approximation 
is  defined  as 


-  flight  phase  - ►j^ —  support  phase 

toe  off  foot  strike 

Figure  1  Human  running  cycle  (Photo  by  Edward  Muybridge) 


OPTIMIZATION  -BASED  PREDICTIVE  DYNAMICS 

Using  optimization  techniques,  the  digital  human’s 
motion  can  be  predicted  as  along  with  the  relative  joint 
torques  and  external  forces.  The  basic  idea  is  to 
determine  joint  angle  profiles  and  torque  profiles  to 
optimize  some  objective  function  (for  example,  metabolic 
energy  consumption  and  joint  torques).  Figure  2  explains 
the  entire  optimization  process.  First,  the  initial  control 
point  values  for  the  joint  angle  profiles  are  given.  Then, 
the  control  points  are  passed  on  to  the  analysis  module. 
The  analysis  module  uses  the  B-spline  module,  DH 
module,  equations  of  motion  module,  and  cost 
function/constraints  module.  Through  this  module,  cost 
function  and  constraints  values  are  obtained.  Using  the 
cost  function  and  constraint  functions  values  and  their 
gradients,  the  optimization  module  checks  whether  or 
not  the  current  values  are  optimum.  The  SNOPT 
program  is  used  for  the  optimization  module.  If  the 
stopping  criteria  are  not  satisfied,  the  B-spline  control 
points  are  updated  by  the  optimization  module.  The 
entire  process  is  repeated  again  until  an  optimum 
solution  is  obtained.  The  B-spline  method,  DH  method, 
and  equations  of  motion  are  discussed  in  the  next 
section. 


q(t)  =  YJNi,p(t)qi  (1) 

i= 0 

where  N,p  is  the  ;th  basis  function  of  degree  p.  The  basis 
function  Nip(t)  is  given  in  the  recursive  form  as 


Nt,p(t,  t)  = 


V  t‘+p  r>  j 


Ni,P- i(Tt)  + 


^  t  —  t^ 

li+p+l  ' 
y^i+p+i  ~  ^i+i  j 


(2) 


and 


Nt,  0(f,t) 


l  (t;  <t<tM) 
0  otherwise 


(3) 


The  relation  between  the  number  of  knots  m+ 1  and  the 
number  of  control  points  n+ 1  is 

m  =  n  +  p  + 1  (4) 


CUBIC  B-SPLINE  CURVES 

The  Cubic  B-spline  curves  which  have  basis  functions  of 
degree  3  in  the  local  interval  t.  <  t  <  tM  is 


q(t)  = 

7=0 


>3 


(5) 


and 


^-3,3  = 


(?i+ 1  h  )(^i+ 1  h-\  )(^'+ 1  h-i ) 


(6a) 


Figure  2  Optimization-based  Predictive  Dynamics  Process 

JOINT  ANGLE  PROFILE  APPROXIMATION 

B-SPLINE  APPROXIMATION 


(ft+ 1 —  ^-2X^+1 —  (-1  Xr,+2  ) 

fX+i  h  XX+2  —  h-\  X^+2  —  h  ) 


(6b) 


fX;+i  —  h- 1  )(^'+ 1  —  h  )(^i+ 2  —  h-\ ) 
|  (t  -  -r,)(r,- +2-r) 

(h+i  —  h  )(^'+ 2  —  h-\  )(^'+ 2  —  h  ) 

,  (?-f,)2(f,+3-0 


The  joint  angles  are  functions  of  time.  These  functions 
can  be  represented  as  a  linear  combination  of  cubic  B- 
spline  basis  functions.  Given  a  knot  vector  t  =  {t0,  tu 


N: 


(t-t,)3 

(h+ 1  —  h  )(tj+ 2  ~  )(f/+3  —  tj ) 


(6d) 


(8) 


To  generate  clamped  B-spline  curves  which  touch  the 
first  and  last  control  points,  multiple  knots  of  multiplicity  4 
are  used  at  the  first  and  the  last  knots. 

KINEMATIC  MODEL  OF  HUMAN  BODY 

DENAVIT-HARTENBERG  METHOD 

Denavit  and  Hartenberg  (1955)  proposed  a  matrix 
transformation  method  to  describe  the  translational  and 
rotational  relationship  systematically  between  adjacent 
links  in  articulated  chain.  This  matrix  transformation 
representation  is  called  the  DH  method.  The 
transformation  matrix  is  a  4x4  homogeneous  matrix. 
This  method  represents  each  link  coordinate  system  in 
terms  of  the  previous  link  coordinate  system.  Then  any 
local  coordinate  system  (including  the  end-effector  of  the 
manipulator  or  serial  chain)  can  be  expressed  in  an 
original  reference  by  the  DH  method.  Basically,  the 
method  represents  a  vector  in  one  coordinate  frame  in 
terms  of  the  other  coordinate  frame.  This  method  has  its 
base  in  the  field  of  robotics  but  can  be  used  for  modeling 
human  kinematics  as  well. 

Consider  articulated  chains,  which  are  depicted  in  Figure 
3. 


Figure  3  Articulated  chains 

Any  point  of  interest  in  the  ith  frame  ;x  can  be 
transferred  to  the  global  reference  frame  °x : 

°x=°i;'x  (7) 

where  'x  is  a  4x1  vector  in  terms  of  the  ith  reference 
frame  and  °T;  is  a  4x4  homogeneous  transformation 
matrix  from  the  *  reference  frame  to  the  global 
reference  frame. 

Here  the  transformation  of  a  vector  to  the  global 
reference  frame  is  simply  multiplication  of  transformation 
matrices,  which  is  given  as: 


°T  =  °T,  'T2---  ,hT  =PJ'T 


The  transformation  matrix  of  this  vector  is  a  4x4  matrix 
that  includes  4  DH  parameters,  which  are  described  in 
Figure  4. 


Figure  4  DH  parameters 

DH  parameters  in  Figure  4  are  defined  as  follows: 

•  9i  is  the  joint  angle  between  the  xM  axis  and  the  x, 
axis  about  the  zM  axis  according  to  the  right-hand 
rule. 

•  dj  is  the  distance  between  the  origin  of  the  /-1th 
coordinate  frame  and  the  intersection  of  the  zM  axis 
with  the  x,  axis  along  the  zM  axis. 

•  a,  is  the  distance  between  the  intersection  of  the  zM 
axis  with  the  x,  axis  and  the  origin  of  the  ith  frame 
along  the  x,  axis.  Or,  the  shortest  distance  between 
the  zM  and  z,  axes. 

•  a,  is  the  angle  between  the  zM  axis  and  the  z,  axis 
about  the  x,  axis  according  to  the  right-hand  rule. 

Then,  the  DH  transformation  matrix  from  ith  frame  to 

/-1th  frame  is  written  as: 


‘T.  = 


cos  6* 
sin  Oi 
0 
0 


-cos or  sin  0j 
cos  or.  cos  6j 
sin  at 
0 


sin  or.  sin  0 
-sin  ar  cos  0t 
cos  a;. 

0 


a.  cos  0j 
aj  sin  6j 
d, 

1 


(9) 


Among  the  four  DH  parameters,  the  rotation  about  the  z 
axis  is  treated  as  the  rotational  DOF  in  the  mechanical 
model,  and  the  other  three  parameters  are  fixed. 
Therefore,  each  transformation  has  one  DOF.  The 
current  Santos™  model  has  55  DOF,  including  6  global 
DOF.  The  global  DOF  are  composed  of  three 
translations  and  three  rotations.  A  full-body  digital 
human  model  is  depicted  in  Figure  5.  Note  that  spine, 
neck,  shoulder  and  hip  joint  have  3  rotational  DOF. 


(11b) 


Elbow,  clavicle,  ankle  and  wrist  joint  have  2  rotational 
DOF.  Knee  and  toe  joint  have  1  rotational  DOF. 


Figure  5  Mechanical  structure  of  Santos  based  on  DH  method 

KINEMATIC  ANALYSIS  OF  HUMAN  BODY 

The  kinematics  analysis  in  the  recursive  form  leads  to  a 
simpler  form  for  the  transformation  matrix  A,  Time 
derivatives  of  the  transformation  matrix  A,  can  be 
obtained  in  the  recursive  form  as: 

A,=AmT,  (10a) 

Bj=A,.=B,._1T,.+A,._1|?H  (10b) 
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DYNAMIC  EQUATIONS  OF  MOTION 
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where  q  is  the  joint  angle  and  T,  is  the  link  transformation 
matrix.  The  gradients  of  transformation  matrices  with 
respect  to  joint  angles,  joint  angle  velocities,  and  joint 
angle  accelerations  are 
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Dynamic  equations  of  motion  are  important  constraints 
in  the  optimization-based  predictive  dynamics  problem  of 
human  running.  The  biggest  challenge  is  the  number  of 
calculations  to  be  performed  because  there  are  many 
matrix  multiplications  and  additions  for  kinematics 
analysis.  Also,  the  optimization  process  can  take  several 
iterations.  These  issues  will  be  discussed  in  this  section. 
Uicker  (1965)  derived  the  standard  formulation  for 
manipulator  dynamics  based  on  Lagrangian  dynamics 
using  DH  4x4  matrix  transformations.  However,  that 
formulation  takes  order  n  calculations.  In  1979,  Waters 
noticed  that  a  simpler  formulation  can  be  derived  that 
takes  only  order  n  calculations.  After  that,  Hollerbach 
(1980)  derived  a  recursive  formulation  from  the  Waters 
formula  that  takes  only  order  n  calculations.  Since  we 
are  solving  an  optimization  problem,  the  number  of 
multiplications  and  additions  that  need  to  be  performed 
are  significant. 


RECURSIVE  LAGRANGE  DYNAMICS  FORMULATION 


(14c) 


The  Lagrange’s  equation  is  given  as 
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where  L  =  K-V  (kinetic  energy  -  potential  energy),  q  is 
the  generalized  coordinate  vector  (joint  angles),  and  r,  is 
joint  torque  vector.  If  any  non-conservative  force  exists, 
it  goes  to  the  left  side  of  Eq.  (12).  When  f  and  h  are 
conservative  external  global  force  and  moment  vectors 
acting  on  the  linkage,  respectively,  Eq.  (12)  can  be 
transformed  to  a  recursive  form  given  as 


D,  =  J,Af  +Tj+1Di+1  (13a) 

E;  =  mi  ‘r;  +  T)+1E;+1  (13b) 

F,  =  V,  +^F,+I  (13c) 

G,=hA+G;+1  (13d) 
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COMPUTATIONAL  CONSIDERATION 

The  number  of  multiplications  and  additions  for  each 
formulation  are  summarized  in  Table  1.  The  order  of 
calculations  for  the  three  formulations  noted  previously 
can  be  observed  in  the  table.  For  a  system  with  small 
DOF,  the  total  computational  time  with  the  three 
formulations  may  not  be  too  different.  However,  for  a 
model  with  a  large  number  of  DOF  (such  as  the  Santos 
model  with  55  DOF),  the  number  of  calculations  can  be 
significantly  different.  This  can  have  a  significant  impact 
on  the  efficiency  of  the  entire  optimization  process.  It  is 
clear  that  the  recursive  formulation  is  the  most  suitable 
for  digital  human  modeling,  and  it  is  used  for  the  running 
problem. 


Table  1  Number  of  multiplication  and  additions 
n:  number  of  transformation  matrices) _ 


Method 

Multiplications 

Additions 

Uicker  (1965) 

32  1/2n4+86  5/12 
n3 

+171  1/4  n2  +  5 
1/3  n  -128 

25  n 4  +  66  1/3  n3 
+129  1/2  n2  +42 
1/3  n  -96 

Waters  (1979) 

106  1/2  ~r?  +  620 

1/2  n-  512 

82  n2  + 514n- 
384 

Hollerbach  (1980) 

830  n -  592 

675  n -  464 

where  J,  is  the  inertia  matrix  for  link  i,  g  is  gravity  vector, 
'r,  is  the  location  of  the  center  of  mass  in  the  ith  local 
frame,  krf  is  the  location  of  the  external  force  acting  in 
the  kth  frame,  z0  =  [0  0  1  0]r,  and  6ik  is  Kronecker  delta. 
The  segment  masses  in  the  mechanical  model  are 
calculated  using  a  mass  distribution  formula  (Chaffin  and 
Andersson,  1984).  The  link  length  and  joint  locations  are 
determined  based  on  high  resolution  3D  scanned  data 
(Eyetronics).  All  the  segments  are  assumed  as  slender 
bars,  and  the  mass  moments  of  inertia  are  calculated 
under  this  assumption. 

The  derivatives  of  equations  of  motion  with  respect  to 
joint  angles,  joint  angle  velocities,  and  joint  angle 
accelerations  are 


Table  2  Number  of  multiplications  and  additions  for 
n= 55 


Method 

Multiplications 

Additions 

Total 

Uicker  (1965) 

312293722 

240195803 

552489525 

Waters  (1979) 

355778 

275936 

631714 

Hollerbach 

(1980) 

45058 

36661 

81719 

VERIFICATION  OF  EQUATIONS  OF  MOTION 

The  equations  of  motion  were  verified  by  the  forward 
dynamics  process  using  a  commercial  general-purpose 
multi-body  dynamics  software  code  (ADAMS).  The 
single  pendulum  problem  was  solved  for  this  process. 
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Figure  6  Single  pendulum 


Figure  6  depicts  the  single  pendulum  model  in  which 
mass  is  0.5  kg,  length  is  0.4  m,  and  it  is  assumed  to  be  a 
slender  bar.  The  equation  of  motion  is  given  by 


In  Figure  8,  point  D  is  ZMP,  which  needs  to  be 
determined.  The  resultant  moment  about  the  ZMP  by 
inertia,  gravity,  and  external  force  (IGF)  is  given  as 


Iq  +  mg—cosq  =  T 


(15) 


MDF  =XDGXmg -XDGXm*G  HC  (16) 


The  initial  position  is  q=  0,  and  the  ADAMS  results  are 
shown  in  Figure  7. 


where  HG  is  rate  of  angular  momentum  about  the 

center  of  mass  of  the  system.  The  resultant  moment 
about  the  point  O  is 


Single  Pendulum 


Figure  7  Join  angle  and  joint  angle  velocity  of  single  pendulum 

Since  it  is  a  free  vibrations  problem,  there  should  not  be 
non-conservative  joint  torque  (r  =  0).  By  using  ADAMS, 
the  joint  torque  is  indeed  obtained  as  zero,  thus  verifying 
the  current  equations  of  motion  formulation. 


M0F  =XOGXmSG-XOGXmiG-^G  (17) 

Then  Eq.  (17)  can  be  written  as 

M„f  =  M®f  - x0D  x R,GF  (18) 

From  the  condition  that  the  tripping  moment  by  the  IGF 
measured  at  the  D  is  zero,  we  have 

nxM'f  =  nxMGGF-nx(x0£)xR,GF) 

=  nxM'CF-(n.R"  )OD  +  (n  •  x0D  )R,GF  (19) 

=  0 


where  n  is  a  unit  vector  that  is  normal  to  ground  plane. 
Then,  the  ZMP  location  is  obtained  as 


_  nxMGGF 
~  n  •  R/gf 


(20) 


STABILITY 


ZERO  YAWING  MOMENT 


ZERO  MOMENT  POINT 

Another  important  consideration  for  the  running  problem 
is  the  dynamic  stability  of  the  motion.  The  most  common 
constraint  to  achieve  stability  for  biped  gait  analysis  is 
the  zero  moment  point  (ZMP)  constraint  in  the  support 
phase.  Zero  moment  point  can  be  derived  by  using  the 
following  steps. 


The  zero  yawing  moment  (ZYM)  constraint  is  usually 
imposed  for  the  upper-body  motion  to  be  compensated 
by  the  lower-body  motion.  From  Eq.  (16),  the  yawing 
moment  about  the  ZMP  D  is  obtained  as 

Ytf!F  =M'!;r  n 

nbody  T  .  "I 

=  I  [x DG.  x  (m,g  -  m,;XGi )  -  Hf.  J-  n  (21) 

i 

nbody  r  -I 

=  I  |x  DG  x  ( mi  g  -  m;x  G ;.)  J  •  11 


where  Hc  is  assumed  to  be  zero. 

FORMULATION 

The  problem  is  to  determine  the  joint  angle  profiles  that 
minimize  an  energy  cost  function.  It  is  assumed  that  the 
running  motion  is  completely  periodic  and  symmetric  and 
that  there  are  two  phases,  support  and  flight.  To  solve 
this  optimization  problem,  a  skeletal  model  of  the  human 
and  the  running  speed  are  needed  as  input.  Through  the 
optimization  process,  joint  angle  profile,  joint  torque 
profile,  and  contact  force  profile  are  obtained  as  output, 
as  shown  in  Figure  9. 


Human 

model 


Running 

speed 


Figure  9  Input  and  output  of  running  problem 

Design  variables  are 
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DV:  q  (22) 

where  q  is  joint  angle  profiles.  The  cost  function  is 

f=\'/'dt  (23) 

which  is  the  proportional  to  the  mechanical  energy.  This 
mechanical  energy  is  a  reasonable  criterion  to  minimize 
(Roussel  et  al.  1998). 

CONSTRAINTS 

Most  constraints  are  motivated  by  the  digital  human 
walking  formulation  (Xiang  et  al.  2007).  The  constraints 
are  listed  as  follows: 


Here,  the  zero  moment  point  is  simply  a  point  where  the 
moments  about  the  x  and  z  axes  due  to  IGF  are  zero. 

Zero  yawing  moment  constraint 

The  yawing  moment  constraint  is  imposed  as 

\y’gf\<yud  (26) 

where  Y™F  is  the  resultant  yawing  moment  about  the  y 
axis  and  T„  is  an  upper  bound  for  it.  In  the  current 
implementation,  TG  is  set  to  zero.  From  Eq.  (21),  the 

zero  yawing  moment  constraint  is  simplified  with  n  =  [0, 
1 , 0] '  as 


1.  Joint  limits 

2.  Ground  penetration 

3.  Foot  location  of  ground  contact  point 

4.  Impact  constraint  (zero  velocity  at  foot  strike) 

5.  ZMP  during  support  phase 

6.  ZYM 

7.  Symmetry  condition 

Current  joint  angle  limits  for  the  body  are  determined 
based  on  Norkin  and  White  (2003). 


nbody 

YDF  =  Z  mi  [(  ~Zzmp)Xi  ~{Xi  ~  xmp  )  \  ]  (27) 

i 

STEP  LENGTH  AND  FLIGHT  TIME 

The  step  length  and  flight  time  were  formulated  as  a 
function  of  running  speed  and  running  frequency, 
respectively  (Bruderlin  and  Calvert,  1996).  The  step 
length  si  is  given  as 


Impact  constraint  (zero  velocity  at  foot  strike) 

As  we  know,  there  is  a  flight  phase  in  human  running.  At 
the  end  of  this  flight  phase,  there  is  impact.  In  this 
impact,  there  is  the  sudden  change  of  joint  angle 
velocities.  Therefore,  this  sudden  change  of  joint  angle 
velocities  results  in  an  impulsive  force  at  the  foot  impact. 
To  handle  the  impact  stage  in  the  current 
implementation,  we  set  the  heel  velocity  to  zero  when 
the  foot  strike  occurs. 

x,.  (t)  —  0,  0  <  t  <  T,  i  e  contact  (24) 

ZMP  constraint 


si  =  0.1394  +  (0.00465  +  level)v^^^^-  (28) 


where  v  is  running  speed  (m/min),  level  is  the  level  of 
expertise  in  running  (-0.001  as  poor  <  level  <  0.001  as 
skilled),  body _height  is  the  height  of  the  human  body.  The 
flight  time  tflight  is  given  as 


tfugh,  =  “0.675x10  3  - (0. 15x10  3  +level)sf 

+  0.542x10 ~V2 


(29a) 


tflighl  =  -8 .925  +  (0. 1 3 1  +  level)  sf  -  0.623  x  1 0-3  sf 

+  0.979xl0‘V3 


(29b) 


To  implement  the  zero  moment  point  constraint  in  the 
current  formulation,  we  consider  the  x-z  plane  as  the 
ground  in  Figure  6.  In  other  words,  the  normal  vector  n  is 
[0,  1,  0]T.  In  this  case,  we  can  simplify  the  ZMP 
calculation  from  Eq.  (20)  as 


where  sf  is  step  frequency  (steps/min,  sf  =  v  I  si).  Eq. 
(29a)  is  used  when  sf  is  0-180  steps/min,  and  Eq.  (29b) 
is  used  when  sf  is  180-230  steps/min. 


RESULTS 


To  evaluate  the  formulation,  models  with  and  without 
arms  were  used.  The  model  without  arms  has  26  DOF  (6 
global  DOF,  7  DOF  for  each  leg,  and  6  DOF  for  spine). 
Figure  10(a)  depicts  joints  in  the  model  without  arms, 
and  Figure  10(b)  depicts  joints  in  the  full-body  model  (55 
DOF). 


Figure  10  (a)  Model  without  arms  (b)  Full-body  model 


The  number  of  control  points  is  taken  as  5  for  each  DOF. 
Thus,  the  total  number  of  design  variables  is  130  for  the 
model  without  arms  and  275  for  full-body  model.  An  Intel 
Pentium  3.46  GHz  CPU  PC  was  used  to  obtain  the 
optimum  solutions. 

Figure  1 1  is  a  snapshot  of  Santos™  running  at  a  speed 
of  2  m/s.  The  step  length  is  0.8  m,  and  the  model  without 
the  arms  was  used  for  this  simulation. 


Figure  11  Snapshot  of  Santos  running  (model  without  arms) 


Figure  12  is  a  snapshot  for  the  case  where  a  backpack 
is  included.  The  model  without  arms  was  used  for  this 
case  as  well.  The  running  speed  was  1.8  m/s  and  step 
length  was  0.6  m.  The  backpack  mass  was  10.20  kg 
(100  N). 


Figure  12  Snapshot  of  Santos  running  with  backpack  (model 
without  arms) 


Figure  13  is  a  snapshot  of  Santos  running  with  the  full- 
body  model.  In  this  simulation,  the  initial  and  end  points 
were  specified  for  the  elbow  as  additional  constraints. 


Figure  13  Snapshot  of  Santos  running  with  arm  motion 

Figure  14  is  a  comparison  of  the  joint  angles  for  the 
spine  between  normal  running  and  running  with 
backpack  (Figures  11  and  12). 


Figure  15  and  Figure  16  are  a  right  knee  joint  angle 
profile  and  a  ground  reaction  force  profile  respectively 
for  the  full-body  model. 


Figure  15  Right  knee  joint  angle  profile 


Figure  16  Ground  reaction  force  on  right  foot 

The  simulation  results  are  compared  those  from  the 
experiments  (Patla  et  al.  1989,  Ounpuu  1994).  The  two 
results  do  not  match  exactly  since  the  dimension  and 
mass  properties  of  the  experimental  subjects  are 
different  from  those  for  the  mechanical  model.  However, 
the  trend  are  similar  to  the  simulation  results.  For  the 
ground  reaction  force,  shape  of  the  impact  moment  is  a 
little  different  from  the  experimental  data.  These  aspects 
need  to  be  investigated  further  to  refine  the  formulation 
of  the  problem. 

CONCLUSION 

The  task  of  digital  human  running  was  formulated  as  an 
optimization  problem.  Using  the  optimization  process,  it 
is  possible  to  predict  dynamic  motion  (joint  angle 
profiles)  as  well  as  the  corresponding  joint  torques.  A 
predictive  dynamics  approach  was  used  where  there 
was  no  need  to  integrate  the  equations  of  motion,  as 
with  the  forward  dynamics  formulation.  B-spline 
interpolation  was  used  for  discretization  along  the  time 
axis,  and  the  Denavit-Hartenberg  method  was  used  for 
kinematics  analysis  of  the  mechanical  system.  For 
dynamic  equilibrium,  the  recursive  Lagrange  method 
was  used  to  reduce  the  order  of  computations.  For 
dynamic  stability,  zero  moment  point  and  zero  yawing 
moment  constraints  were  used.  To  formulate  the  impact 


stage,  the  zero  velocity  at  foot  strike  was  used.  The 
mechanical  structure  of  Santos™  was  developed  with  (1) 
a  model  without  arms,  which  had  26  DOF,  and  (2)  a  full- 
body  model,  which  had  55  DOF.  As  a  separate  case,  an 
external  load  was  applied  as  a  backpack.  With  the  full- 
body  model,  we  could  observe  the  upper-body  motion, 
especially  the  arm  motion.  The  step  length  and  flight 
time  were  given  as  a  functions  of  running  speed  and 
running  frequency,  respectively.  A  more  detailed 
validation  of  the  formulation  for  the  running  problem  is  in 
progress  that  will  be  reported  later. 
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